One-Dimensional Climate Models: The Basics of Brown Dwarfs¶
In this tutorial you will learn the very basics of running 1D climate runs. For a more in depth look at the climate code check out Mukherjee et al. 2022 (note this should also be cited if using this code/tutorial).
What you should already be familiar with:
What you will need to download to use this tutorial:
Download 1460 PT, 196 wno Correlated-K Tables from Roxana Lupu to be used by the climate code for opacity
Download the sonora bobcat cloud free
structures_file so that you can validate your model run
Note: the two files above are dependent on metallicity and C/O. For this tutorial we will stick to solar M/H and solar C/O, but note that you can change that by picking the right C-K file in the opannection step
[22]:
import warnings
warnings.filterwarnings('ignore')
import picaso.justdoit as jdi
import picaso.justplotit as jpi
jpi.output_notebook()
import astropy.units as u
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
[2]:
#1 ck tables from roxana
mh = '+000'#'+0.0' #log metallicity
CtoO = '100'#'1.0' # CtoO ratio
ck_db = f'/Users/nbatalh1/Documents/data/kcoeff_2020/sonora_2020_feh{mh}_co_{CtoO}.data.196'
#ck_db = f'/Users/sagnickmukherjee/Documents/GitHub/Disequilibrium-Picaso/reference/climate_INPUTS/sonora_2020_feh+1.0_co_1.0.data.196'
#sonora bobcat cloud free structures file
sonora_profile_db = '/Users/nbatalh1/Documents/data/sonora_bobcat/structures_m+0.0/'
What does a climate model solve for?¶
1D Radiative-Convective Equilibrium Models solve for atmospheric structures of brown dwarfs and exoplanets, which includes:
The Temperature Structure (T(P) profile)
The Chemical Structure
Energy Transport in the atmosphere
But these physical components are not independent of each other. For example, the chemistry is dependent on the T(P) profile, the radiative transfer is dependent on clouds and the chemistry and so on.
PICASO tries to find the atmospheric state of your object by taking care of all of these processes and their interconnections self-consistently and iteratively. Therefore, you will find that the climate portion of PICASO is slower than running a single forward model evaluation.
Starting up the Run¶
You will notice that starting a run is nearly identical as running a spectrum. However, how we will add climate=True to our inputs flag. We will also specify browndwarf in this case, which will turn off the irradiation the object is receiving.
New Parameter: Effective Temperature. This excerpt from Modeling Exoplanetary Atmospheres (Fortney et al) provides a thorough description and more reading, if you are interested.
If the effective temperature, \(T_{eff}\), is defined as the temperature of a blackbody of the same radius that would emit the equivalent flux as the real planet, \(T_{eff}\) and \(T_{eq}\) can be simply related. This relation requires the inclusion of a third temperature, \(T_{int}\), the intrinsic effective temperature, that describes the flux from the planet’s interior. These temperatures are related by:”
\(T_{eff} = T_{int} + T_{eq}\)
We then recover our limiting cases: if a planet is self-luminous (like a young giant planet) and far from its parent star, \(T_{eff} \approx T_{int}\); for most rocky planets, or any planets under extreme stellar irradiation, \(T_{eff} \approx T_{eq}\).
[3]:
cl_run = jdi.inputs(calculation="browndwarf", climate = True) # start a calculation
#note you need to put the climate keyword to be True in order to do so
# now you need to add these parameters to your calculation
teff= 1000 # Effective Temperature of your Brown Dwarf in K
grav = 1000 # Gravity of your brown dwarf in m/s/s
cl_run.gravity(gravity=grav, gravity_unit=u.Unit('m/(s**2)')) # input gravity
cl_run.effective_temp(teff) # input effective temperature
Let’s now grab our gaseous opacities, whose path we have already defined above. Again, this code uses a correlated-k approach for accurately capturing opacities (see section 2.1.4; Mukerjee et al 2022).
[4]:
# Notice The keyword ck is set to True because you want to use the correlated-k opacities for your calculation
# and not the line by line opacities
opacity_ck = jdi.opannection(ck_db=ck_db) # grab your opacities
Initial T(P) Guess¶
Every calculation requires an initial guess of the pressure temperature profile. The code will iterate from there to find the correct solution. A few tips:
We recommend using typically 51-91 atmospheric pressure levels. Too many pressure layers increases the computational time required for convergence. Too little layers makes the atmospheric grid too coarse for an accurate calculation.
Start with a guess that is close to your expected solution. We will show an example using an isothermal P(T) profile below so you can see the iterative process. Later though, we recommend leveraging pre-computed grids (e.g. Sonora) as a starting guess for Brown Dwarfs.
[5]:
nlevel = 91 # number of plane-parallel levels in your code
#Lets set the max and min at 1e-4 bars and 500 bars
Pmin = 1e-4 #bars
Pmax = 500 #bars
pressure=np.logspace(np.log10(Pmin),np.log10(Pmax),nlevel) # set your pressure grid
temp_guess = np.zeros(shape=(nlevel)) + 500 # K , isothermal atmosphere guess
Initial Convective Zone Guess¶
You also need to have a crude guess of the convective zone of your atmosphere. Generally the deeper atmosphere is always convective. Again a good guess is always the published SONORA grid of models for this. But lets assume that the bottom 7 levels of the atmosphere is convective.
New Parameters:
nofczns: Number of convective zones. Though the code has functionality to solve for more than one. In this basic tutorial, let’s stick to 1 for now.rfacv: (See Mukherjee et al Eqn. 20r_st) https://arxiv.org/pdf/2208.07836.pdfnstr_upper: this defines the top most level of your guessed convective zone. If you don’t have a clue where your convective zone might end be choose a number that is \(\sim\)nlevel-5 (a few pressure levels away from the very bottom of your grid)
Non-zero values of rst (aka “rfacv” legacy terminology) is only relevant when the external irradiation on the atmosphere is non-zero. In the scenario when a user is computing a planet-wide average T(P) profile, the stellar irradiation is contributing to 50% (one hemisphere) of the planet and as a result rst = 0.5. If instead the goal is to compute a night-side average atmospheric state, rst is set to be 0. On the other extreme, to compute the day-side atmospheric state of a tidally locked planet rst should be set at 1.
[6]:
nofczns = 1 # number of convective zones initially. Let's not play with this for now.
nstr_upper = 83 # top most level of guessed convective zone
nstr_deep = nlevel -2 # this is always the case. Dont change this
nstr = np.array([0,nstr_upper,nstr_deep,0,0,0]) # initial guess of convective zones
# Here are some other parameters needed for the code.
rfacv = 0.0 #we are focused on a brown dwarf so let's keep this as is
Now we would use the inputs_climate function to input everything together to our cl_run we started.
[7]:
cl_run.inputs_climate(temp_guess= temp_guess, pressure= pressure,
nstr = nstr, nofczns = nofczns , rfacv = rfacv)
Run the Climate Code¶
The actual climate code can be run with the cl_run.run command. The save_all_profiles is set to True to save the T(P) profile at all steps. The code will now iterate from your guess to reach the correct atmospheric solution for your brown dwarf of interest.
[8]:
out = cl_run.climate(opacity_ck, save_all_profiles=True,with_spec=True)
Iteration number 0 , min , max temp 499.98326466392086 837.9263214798718 , flux balance -0.9374419003155262
Iteration number 1 , min , max temp 499.90282151307053 1128.933155611243 , flux balance -0.936408261722275
Iteration number 2 , min , max temp 499.0739010829525 1989.8539978094143 , flux balance -0.7652890988601658
Iteration number 3 , min , max temp 488.3467473257478 2051.0908072121306 , flux balance -0.6761972266034364
Iteration number 4 , min , max temp 461.13324815007417 2161.3900258729986 , flux balance -0.4562694693726309
Iteration number 5 , min , max temp 365.3039004488172 2351.913366019811 , flux balance 0.16976424728623773
Iteration number 6 , min , max temp 305.1084929609062 2318.798344605897 , flux balance 0.017634000036814805
Iteration number 7 , min , max temp 273.94640826204846 2314.9720483881833 , flux balance 0.0009617689167044826
Iteration number 8 , min , max temp 264.8032004487643 2314.7331025689386 , flux balance 5.218402169497495e-05
Iteration number 9 , min , max temp 264.0688523260907 2314.719600921175 , flux balance 2.9239156762267505e-06
In t_start: Converged Solution in iterations 9
Big iteration is 264.0688523260907 0
Iteration number 0 , min , max temp 266.65590959585927 2432.127389376785 , flux balance -0.7970694548689117
Iteration number 1 , min , max temp 276.38631835941004 2774.853389806402 , flux balance -0.6990886542916896
Iteration number 2 , min , max temp 291.585259943468 3068.7646336149924 , flux balance -0.5300844485358708
Iteration number 3 , min , max temp 342.95658573476436 3731.5842564085374 , flux balance 0.34440717195673387
Iteration number 4 , min , max temp 333.09150636381 3569.1500194802065 , flux balance 0.03610415718845372
Iteration number 5 , min , max temp 332.1938282718999 3551.0906913649633 , flux balance 0.0005505219829357028
Iteration number 6 , min , max temp 332.1856522339841 3550.8884923979317 , flux balance 2.7336245585222263e-07
In t_start: Converged Solution in iterations 6
Big iteration is 332.1856522339841 1
Iteration number 0 , min , max temp 332.7405410985925 3759.4200541604387 , flux balance -0.179254805861405
Iteration number 1 , min , max temp 339.48470946215497 5199.9 , flux balance -0.12087123748731982
Iteration number 2 , min , max temp 339.8032137583409 5199.9 , flux balance -0.11800464398944968
Iteration number 3 , min , max temp 340.1034378619028 5199.9 , flux balance -0.1152941973750462
In t_start: Converged Solution in iterations 3
Big iteration is 340.1034378619028 2
Iteration number 0 , min , max temp 340.1079050774916 5199.9 , flux balance -0.12324068523678665
In t_start: Converged Solution in iterations 0
Profile converged
Iteration number 0 , min , max temp 340.1123697488552 5199.9 , flux balance -0.12320712699273484
In t_start: Converged Solution in iterations 0
Big iteration is 340.1123697488552 0
Iteration number 0 , min , max temp 340.11683181904954 5199.9 , flux balance -0.1231735873891928
In t_start: Converged Solution in iterations 0
Profile converged
Move up two levels
Iteration number 0 , min , max temp 341.9660328163306 5013.6051527766795 , flux balance -0.10726949162612438
Iteration number 1 , min , max temp 354.30268332901954 5199.9 , flux balance 0.006182701340163212
Iteration number 2 , min , max temp 353.6061965104687 5199.9 , flux balance 2.520625619959532e-05
Iteration number 3 , min , max temp 353.60363903812646 5199.9 , flux balance 6.715077108058855e-09
In t_start: Converged Solution in iterations 3
Big iteration is 353.60363903812646 0
Iteration number 0 , min , max temp 353.65461173450444 5199.9 , flux balance -0.017681682545377667
Iteration number 1 , min , max temp 353.79690682937814 5199.9 , flux balance -0.016874443201559366
Iteration number 2 , min , max temp 353.90654681622794 5199.9 , flux balance -0.016251735980188527
Iteration number 3 , min , max temp 353.99522308686113 5199.9 , flux balance -0.015747651135034485
Iteration number 4 , min , max temp 354.00264098871367 5199.9 , flux balance -0.0157054725168527
In t_start: Converged Solution in iterations 4
Big iteration is 354.00264098871367 1
Iteration number 0 , min , max temp 354.0068257833241 5199.9 , flux balance -0.016263467988756213
In t_start: Converged Solution in iterations 0
Profile converged
Move up two levels
Iteration number 0 , min , max temp 355.6751719778109 4812.132539130671 , flux balance -0.006533293437661223
Iteration number 1 , min , max temp 356.7991638335737 4934.410979595851 , flux balance 5.3015617592315145e-05
Iteration number 2 , min , max temp 356.79316381155945 4931.255778229213 , flux balance 2.088976667215314e-08
In t_start: Converged Solution in iterations 2
Big iteration is 356.79316381155945 0
Iteration number 0 , min , max temp 357.2502527549278 5125.317089426997 , flux balance 5.4451235506153185e-05
Iteration number 1 , min , max temp 357.24920692007095 5118.009062814578 , flux balance 4.949314728709674e-08
In t_start: Converged Solution in iterations 1
Big iteration is 357.24920692007095 1
Iteration number 0 , min , max temp 357.3250039864804 5199.9 , flux balance 1.2052355146574003e-05
Iteration number 1 , min , max temp 357.3249665337781 5199.9 , flux balance 4.770254500824475e-09
In t_start: Converged Solution in iterations 1
Big iteration is 357.3249665337781 2
Iteration number 0 , min , max temp 357.3429191665995 5199.9 , flux balance 2.230408274636585e-06
Iteration number 1 , min , max temp 357.34291618913136 5199.9 , flux balance 4.827892634990887e-10
In t_start: Converged Solution in iterations 1
Big iteration is 357.34291618913136 3
Iteration number 0 , min , max temp 357.3483687517044 5199.9 , flux balance 3.0147324869194343e-07
Iteration number 1 , min , max temp 357.348368606285 5199.9 , flux balance 5.155788103969342e-11
In t_start: Converged Solution in iterations 1
Big iteration is 357.348368606285 4
Not converged
Move up two levels
Iteration number 0 , min , max temp 357.3500686276559 4072.364990266688 , flux balance 3.4297795606697915e-08
In t_start: Converged Solution in iterations 0
Big iteration is 357.3500686276559 0
Iteration number 0 , min , max temp 357.3507633432517 4073.0074447077304 , flux balance 7.487024525653784e-10
In t_start: Converged Solution in iterations 0
Profile converged
Move up two levels
Iteration number 0 , min , max temp 357.41985712242104 3527.0855534636707 , flux balance 1.1961114355037407e-05
In t_start: Converged Solution in iterations 0
Big iteration is 357.41985712242104 0
Iteration number 0 , min , max temp 357.44804239229796 3529.969168241415 , flux balance 9.96211577753417e-07
In t_start: Converged Solution in iterations 0
Profile converged
Iteration number 0 , min , max temp 357.62940680268235 3473.0083624064778 , flux balance 6.257157102743745e-05
Iteration number 1 , min , max temp 357.6292589504265 3472.8552394653148 , flux balance 7.947525569256054e-08
In t_start: Converged Solution in iterations 1
Big iteration is 357.6292589504265 0
Iteration number 0 , min , max temp 357.7270229820692 3476.5066458097835 , flux balance 3.3446028886451727e-06
In t_start: Converged Solution in iterations 0
Profile converged
Iteration number 0 , min , max temp 357.86095978055914 3463.4849206372023 , flux balance 1.0082507971320976e-05
In t_start: Converged Solution in iterations 0
Big iteration is 357.86095978055914 0
Iteration number 0 , min , max temp 357.9430737911975 3464.980597100013 , flux balance -1.9523117208982595e-06
In t_start: Converged Solution in iterations 0
Profile converged
final [ 0 73 89 0 0 0]
Iteration number 0 , min , max temp 357.97336872109037 3465.3460498763943 , flux balance -7.562840261377672e-07
In t_start: Converged Solution in iterations 0
Big iteration is 357.97336872109037 0
Iteration number 0 , min , max temp 357.98220214344207 3465.4356049224984 , flux balance -2.0350508495295448e-07
In t_start: Converged Solution in iterations 0
Profile converged
YAY ! ENDING WITH CONVERGENCE
Benchmark with Sonora Bobcat¶
[9]:
pressure_bobcat,temp_bobcat = np.loadtxt(jdi.os.path.join(
sonora_profile_db,f"t{teff}g{grav}nc_m0.0.dat"),
usecols=[1,2],unpack=True, skiprows = 1)
plt.figure(figsize=(10,10))
plt.ylabel("Pressure [Bars]", fontsize=25)
plt.xlabel('Temperature [K]', fontsize=25)
plt.ylim(500,1e-4)
plt.xlim(200,3000)
plt.semilogy(out['temperature'],out['pressure'],color="r",linewidth=3,label="Our Run")
plt.semilogy(temp_bobcat,pressure_bobcat,color="k",linestyle="--",linewidth=3,label="Sonora Bobcat")
plt.minorticks_on()
plt.tick_params(axis='both',which='major',length =30, width=2,direction='in',labelsize=23)
plt.tick_params(axis='both',which='minor',length =10, width=2,direction='in',labelsize=23)
plt.legend(fontsize=15)
plt.title(r"T$_{\rm eff}$= 1000 K, log(g)=5.0",fontsize=25)
[9]:
Text(0.5, 1.0, 'T$_{\\rm eff}$= 1000 K, log(g)=5.0')
Climate Plots and Animations¶
Animate Convergence¶
We can also try to see how our initial guess of an isothermal atmosphere was changed by the code to reach the converged solution
[10]:
ani = jpi.animate_convergence(out, cl_run, opacity_ck,
molecules=['H2O','CH4','CO','NH3'])
[11]:
ani
[11]:
Brightness Temperature¶
Checking the brightness temperature serves many useful purposes:
Intuition building. Allows you to see what corresponding temperature are you sensitive to at each wavelength
Note that this temperature doesn’t need to be the physical temperature of your atmosphere but if you can find the physical converged atmospheric temperature closest to this brightness temperature you can also get an idea of the atmospheric pressure from where the flux you are seeing is originating from.
Determining if your choice in bottom boundary pressure grid was correct.
If your brightness temperature is such that you bottom out at the temperature corresponding to the highest pressure, you have not extended your grid to high enough pressures.
Brightness Temperature Equation:
\(T_{\rm bright}=\dfrac{a}{{\lambda}log\left(\dfrac{{b}}{F(\lambda){\lambda}^5}+1\right)}\)
where a = 1.43877735x\(10^{-2}\) m.K and b = 11.91042952x\(10^{-17}\) m\(^4\)kg/s\(^3\)
Let’s calculate the brightness temperature of our current run and check if our pressure grid was okay.
[12]:
brightness_temp, figure= jpi.brightness_temperature(out['spectrum_output'])
In the above plot you can see that your brightness temperature is nicely bound between the minimum and maximum temperature. Your run is good and your choice of pressure grid is also great. Well done team!
Selecting an Adequate Pressure Grid¶
For understanding and intuition building, let’s do out a run where we purposely choose an incomplete pressure grid. Let’s do the same run by the max pressure set at only 3 bars instead of 500 bars.
[13]:
cl_bad_pres = jdi.inputs(calculation="brown", climate = True)
cl_bad_pres.gravity(gravity=grav, gravity_unit=u.Unit('m/(s**2)')) # input gravity
cl_bad_pres.effective_temp(teff)
nlevel = 91 # number of plane-parallel levels in your code
Pmin = 1e-4 #bars
Pmax = 3 #bars
pressure=np.logspace(np.log10(Pmin),np.log10(Pmax),nlevel) # set your pressure grid
temp_guess = np.zeros(shape=(nlevel)) + 500 # K , isothermal atmosphere guess
nofczns = 1 # number of convective zones initially. Let's not play with this for now.
nstr_upper = 83 # top most level of guessed convective zone
nstr_deep = nlevel -2 # this is always the case. Dont change this
nstr = np.array([0,nstr_upper,nstr_deep,0,0,0]) # initial guess of convective zones
rfacv = 0.0 #we are focused on a brown dwarf so let's keep this as is
cl_bad_pres.inputs_climate(temp_guess= temp_guess, pressure= pressure,
nstr = nstr, nofczns = nofczns , rfacv = rfacv)
out_bad_pres = cl_bad_pres.climate(opacity_ck, save_all_profiles=True,with_spec=True)
Iteration number 0 , min , max temp 496.99226189262305 710.0029306935495 , flux balance -0.8582873586402681
Iteration number 1 , min , max temp 489.494108418519 838.82866668617 , flux balance -0.7427284317668132
Iteration number 2 , min , max temp 482.47475243759106 897.6660760486611 , flux balance -0.6569322383054985
Iteration number 3 , min , max temp 475.4802992080763 936.6610851224333 , flux balance -0.5821469819473322
Iteration number 4 , min , max temp 456.14469172667305 1015.477361029247 , flux balance -0.3653375110906861
Iteration number 5 , min , max temp 398.51105697606516 1117.800883671986 , flux balance 0.14217011056196238
Iteration number 6 , min , max temp 379.1439751136505 1097.639719211233 , flux balance 0.012579232094419022
Iteration number 7 , min , max temp 376.83697980870176 1095.5955242856853 , flux balance 0.00040236012708549956
Iteration number 8 , min , max temp 376.79678083742874 1095.528241574759 , flux balance 1.0464435533938503e-05
In t_start: Converged Solution in iterations 8
Big iteration is 376.79678083742874 0
Iteration number 0 , min , max temp 365.0526662589471 1153.1085250124231 , flux balance 0.034608421931663795
Iteration number 1 , min , max temp 363.8332512844389 1147.0994151287343 , flux balance 0.0014481375714295141
Iteration number 2 , min , max temp 363.78202018184817 1146.8435897729275 , flux balance 3.909661410062007e-05
In t_start: Converged Solution in iterations 2
Big iteration is 363.78202018184817 1
Iteration number 0 , min , max temp 362.93146860391477 1152.4225731682475 , flux balance -0.00037357192013021633
In t_start: Converged Solution in iterations 0
Profile converged
Iteration number 0 , min , max temp 362.9457575575597 1152.8137241650763 , flux balance -5.564146071484752e-05
In t_start: Converged Solution in iterations 0
Big iteration is 362.9457575575597 0
Iteration number 0 , min , max temp 362.9586968663064 1152.8440842883838 , flux balance -4.479567321641991e-06
In t_start: Converged Solution in iterations 0
Profile converged
final [ 0 83 89 0 0 0]
We are already at a root, tolf , test = 5e-05 , 1.576146323072245e-05
Big iteration is 362.9586968663064 0
We are already at a root, tolf , test = 5e-05 , 1.576146323072245e-05
Profile converged
YAY ! ENDING WITH CONVERGENCE
Lets plot the profile from our new run and check it against our old run.
[14]:
plt.figure(figsize=(10,10))
plt.ylabel("Pressure [Bars]", fontsize=25)
plt.xlabel('Temperature [K]', fontsize=25)
plt.ylim(30,1e-4)
plt.xlim(200,1200)
plt.semilogy(out['temperature'],out['pressure'],color="r",linewidth=3,label="Good Run")
plt.semilogy(out_bad_pres['temperature'],out_bad_pres['pressure'],color="b",linewidth=3,label="Bad Pressure Run")
plt.semilogy(temp_bobcat,pressure_bobcat,color="k",linestyle="--",linewidth=3,label="Sonora Bobcat")
plt.minorticks_on()
plt.tick_params(axis='both',which='major',length =30, width=2,direction='in',labelsize=23)
plt.tick_params(axis='both',which='minor',length =10, width=2,direction='in',labelsize=23)
plt.legend(fontsize=15)
plt.title(r"T$_{\rm eff}$= 1000 K, log(g)=5.0",fontsize=25)
[14]:
Text(0.5, 1.0, 'T$_{\\rm eff}$= 1000 K, log(g)=5.0')
This new profile is slightly off from our run and also the sonora bobcat run. Lets look at its brightness temperature as a function of wavelength and check if it matches well with our previous run.
[15]:
brightness_temp_bad, figure_bad= jpi.brightness_temperature(
out_bad_pres['spectrum_output'])
See how the brightness temperature from this new run is different from our previous succesful run. The brightness temperatures infact goes over the maximum temperature achieved by the model. Therefore the pressure grid used for this run is incorrect because one can look through the atmosphere to the bottom of the grid at most wavelengths which is not good and the resultant “converged” T(P) profile is also visibly inaccurate as a result as well.
Post Process High Resolution Spectrum¶
We can quickly do this by resetting the opannection to not use the ck database and use the ptchem_df DataFrame as input for the atmosphere
This is also the point where you could post-process clouds using virga or a box model as seen in these tutorials here: 1. Adding clouds with virga 2. Adding box model clouds
[23]:
opa_mon = jdi.opannection()
hi_res = jdi.inputs(calculation="browndwarf") # start a calculation
teff= 1000 # Effective Temperature of your Brown Dwarf in K
grav = 1000 # Gravity of your brown dwarf in m/s/s
hi_res.gravity(gravity=grav, gravity_unit=u.Unit('m/(s**2)')) # input gravity
hi_res.atmosphere(df=out['ptchem_df'])
df_spec = hi_res.spectrum(opa_mon, calculation='thermal')
w,f = jdi.mean_regrid(df_spec['wavenumber'],df_spec['thermal'], R=100)
[25]:
jpi.show(jpi.spectrum(w,f,x_axis_type='log',y_axis_type='log'))
[ ]: